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We study numerically the critical region and the disordered phase of the random transverse-field 
Ising chain. By using a mapping of Lieb, Schultz and Mattis to non-interacting fermions, we can 
obtain a numerically exact solution for rather large system sizes, L < 128. Our results confirm the 
striking predictions of earlier analytical work and, in addition, give new results for some probability 
distributions and scaling functions. 
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I. INTRODUCTION 

It has recently become clear that quantum phase 
transitions^ in disordered systems are rather different 
from phase transitions-,driven by thermal fluctuations. 
In particular, Griffithsel showed that the free energy is 
a non-analytic function of the magnetic field in part of 
the disordered phase because of rare regions, which are 
more strongly correlated than the average and which are 
locally ordered. However, in a classical system., this effect 
is very weak, all the derivatives being finitea. By con- 
trast, in a quantum system at zero temperature, these 
effects are much more pronounced. 

One model where these effects can be worked out in de- 
tail, and where rare, strongly coupled regions dominate 
not only the disordered phase but also the critical re- 
gion, is the one-dimensional random transverse-field Ising 
chain with Hamiltonian 



"H = — S JiO^o 



i+l 



L 



(i) 



Here the {erf} arc Pauli spin matrices, and the interac- 
tions Ji and transverse fields hi are both independent 
random variables, with distributions 7r(J) and p(h) re- 
spectively. The lattice size is L, which we take to be even, 
and periodic boundary conditions are imposed. The 
ground state of this model is closely related to the finite- 
temperature behavior of a two-dimensional classical Ising 
model with disorder perfectly correlated along one direc- 
tion, which was first studied by McCoy and Will Sub- 
sequently, the quantum model, Eq. ([I]), was studied by 
Shankar and MurphyEl, and recently, in great detail, by 
FisherH. From a real space renormalization group analy- 
sis, which becomes exact on large scales, Fisher obtained 
many new results and considerable physical insight. The 
purpose of the present study is to investigate the,model in 
Eq. numerically, using a powerful techniquaJ which is 
special to one-dimensional systems n to verify the surpris- 
ing predictions of the earlier worHaij and to determine 



certain distributions and scaling functions which have not 
yet been calculated analytically. 

In one-dimension one can perform a gauge transforma- 
tion to make all the Ji and hi positive. Unless otherwise 
stated, the numerical work used the following rectangular 
distribution: 



7T(J) = 



P (h) 



1 for < J < 1 
otherwise 

h^ 1 for < h < h Q 
otherwise. 



(2) 



The model is therefore characterized by a single control 
parameter, h$. As discussed in section II, the critical 
point is at ho — 1 (so the distributions of h and J are 
then the same) and the deviation from criticality is con- 
veniently measured by the parameter S in Eq. (|^), where, 
for the distribution in Eq. (0), 



S = — In hn 
2 



(3) 



Section II discusses the analytical results obtained pre- 
viously, ancLsection Ill-reviews the work of Lieb, Schultz 
and MattisQ, KatsuraO and Pfeutyfl which relates the 
Hamiltonian to free fermions, and also explains how this 
technique can be implemented numerically for the ran- 
dom case. In section IV the numerical results for the 
distribution of the energy gap are shown, while section V 
discusses results for the correlation functions. Results for 
the local susceptibility on smaller sizes, obtained by the 
Lanczos method, are discussed in section VI, while data 
for the q = structure factor, which could be measured 
in a scattering experiment, are considered in section VII. 
Finally, in section VIII, we summarize our conclusions 
and discuss the possible relevance of the results to mod- 
els in higher dimensions. 



II. ANALYTICAL RESULTS 

In this section we summarize the results obtained ear- 
lier by McCoy and WrO, Shankar and MurphyEl and par- 
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ticularly by Fisherfl. Defining 

A h = [lnh] a 
A J = [lnJ] f 



(4) 



where [. . .] av denote an average over disorder, the critical 
point occurs when 



A,, = A.. 



/ ■ 



(•5) 



Clearly this is satisfied if the distributions of bonds and 
fields are equal, and the criticality of the model then fol- 
lows from duality. A convenient measure of the deviation 
from criticality is given by 



6 = 



A h - Aj 



[(ln/i) 2 



Al 



[(In J) 2 



A 2 j 



(6) 



At a quantum critical point one needs to consider the 
dynamical critical exponent, z, even when determining 
static critical phenomena, because statics and dynamics 
are coupled. The relation between a characteristic length 
scale I and the corresponding time scale r is then r ~ V . 
For the present model one has, at the critical point, 



z = oo (5 = 0), 



(7) 



or, more precisely, the time scale varies as the exponential 
of the square root of the corresponding length scale. In 
addition, the distribution of local relaxation times is pre- 
dicted to be very broad. One of the goals of the present 
work is to determine the form of the distribution of a 
related quantity, the gap to the first excited state. 

Moving into the disordered phase, there is still a very 
broad distribution of relaxation times because of Grif- 
fiths singularities, and one can still, as a result, define a 
dynamical exponent but this now varies with 5, diverging 
as 



1 

26 



+ C + 0(5) , 



(8) 



for 8 — > 0, where C is a non-universal constant. Mov- 
ing further away from the critical point, if one reaches a 
situation where all the fields are bigger than all the inter- 
actions, then Griffiths singularities no longer occur and 
the distribution of relaxation times becomes narrow. De- 
noting the value of 8 where this happens by 8q, Griffiths 
singularities occur in that part of the disordered phase 
whereEd 

< 5 < 8 G . (9) 
Approaching the end of the Griffiths phase, one has 



lim z = 



(10) 



For the distribution used in the numerical calculations, 
Eq. (||), Sg — oo so Griffiths singularities occur through- 
out the disordered phase. In the disordered phase, the 



magnetization in the z-direction has a singular piece if a 
uniform field, H, coupling to er z , is added, namely 



'sing 



\H\ 



(11) 



so the linear susceptibility diverges over part of the dis- 
ordered phase, a result first found by McCoy and Wvu. 

Next we turn to predictions for the correlation func- 
tions 



* 3 ' 



(12) 



Again there are very big fluctuations, and, as a result, the 
average and typical correlations behave quite differently. 
The average correlation function, 



1 L 

Cav(r) = ? X>?<H-> 



varies as a power of r at criticality, 



C a v(r) 



1 



(6 = 0), 



where 



1 + V5 



= 1.61804. 



(13) 



(14) 



(15) 



is the golden mean, so the power in Eq. (fL4|) is approxi- 
mately 0.38. Away from criticality, C av (r) decays expo- 
nentially at a rate given by the true correlation length, 
£, where 



with 



v = 2 . 



(16) 



(17) 



The amplitude of the correlation length, ly, is also known 
and given by 



lv = 



(18) 



var [h] + var [J] 
For the distribution in Eq. (0) one has 

ly = 1. (19) 
Scaling theory predicts that 

C a vM) 



C av (r; 6 = 0) 



C av (r/d) 



(20) 



where C av is a universal scaling function and v is the 
true correlation function exponent, predictecfl to equal 
2. Fishem has calculated the asymptotic form of the 
scaling function in Eq. pTj) for r ^> £ and finds 



2 



CM = D- 



-x— 4.055a; 



1/3 



r 0.451 



{X > 1) , 



(21) 



where D is an unknown constant, and 0.451 is the nu- 
merical value of 5/6 — (2 — 0). 

The average correlation function is, however, domi- 
nated by rare pairs of spins which have a correlation func- 
tion of order unity, much larger than the typical value, so 
it is necessary to consider the distribution of In C (r) to 
get an idea of the typical behavior. At the critical point 



In C(r) 



(,5 = 0) 



(22) 



with the coefficient in Eq. (|22j) having a distribution 
which is independent of r. A goal of the present study 
is to investigate this distribution numerically. In the dis- 
ordered phase, — lnC(r) oc r with a coefficient which is 
self-averaging for r — > oo, i.e. 



lnC(r) w r/f 



(23) 



for large r, where the typical correlation length, £, has 
the behavior 



with 



v = 1 



(24) 



(25) 



The scaling equation corresponding to Eq. (|20j) but for 
the average of the log of the correlation function is 



In 



C(r;S)] av 
C(r;<5 = 0) 



lnCi„ p (r/f) 



(26) 



where Ct yp is a universal scaling function. From Eqs. (|22| ) 
and (J23j) one has, for r>(, 



(27) 



For correlations of quantities such as the energy, which 
are local in the fermion operators, see Eqs. (|||) and (p4| ) 
of the next section, Shankar and MurphyO obtained more 
detailed information. They calculated not only the ex- 
ponent for the typical correlation length in Eq. ( p4| ) but 
also the amplitude, finding 



t 1 = [ln/i] av -[lnJ] s 



(28) 



exactly. For r>( the mean of In C en (r) is defined to be 
— r/£ so 



[lnC e „(r)] av « -{[ln/i] av - [In J]] av }r 



(29) 



in this_limit. The variance of the distribution is also 
knowrO for r»^: 



var [In C en (r)] w {var [ln/i] + var [In J]}', 



(30) 



Note that the standard deviation of In C en (r) is propor- 
tional to r 1 / 2 whereas the mean is proportional to r, so 
In C en (rlJpecomes self-averaging for r 3> £. 

FisheiO has suggested that Eqs. (p8|)-(|30|) might also 
be true asymptotically for quantities such as o~ z which 
are not local in fermion operators. If this is true, then, 
for the distribution in Eq. (S), we have 



[lnC(r;6)] av ^-2dr 
var [lnC(r;S)] w 2r 



(31) 
(32) 



for r ^> ^. 

Note that an important feature of these results is that 
the true correlation length (which describes the average 
correlation function) has a different exponent from that 
of the typical correlation length. 



III. MAPPING TO FREE FERMIONS 

The numerical calculations are enormously simplified 
by relating the model in Eq. (^) to non-interacting 
fermions. This technique was first developed for some 
related quantum spin chain problems in a beautiful pa- 
per by Lieb, Schultz and MattisQ and then applied to the 
non-random transverse field Ising chain by Katsuraa and 
PfeutyH. 

The starting point is the Jordan- Wigner transforma- 
tion, which relates the spin operators to fermion creation 
and annihilation operators, cj and Cj, by the following 
transformation: 

rr^ _l_ n ■ 



i(a\ - a,) 
1 - 2aUi = 1 



2c] Ci 



(33) 



where 



cj exp 



exp 



etc 



3=1 
i-1 

J'=l 



(34) 



This works because the Pauli spin matrices anti-commutc 
on the same site but commute on different sites. The 
"string operator" in the exponentials in Eq. ([li]) is just 
what is needed to insert an extra minus sign, converting 
a commutator to an anti-commutator for different sites. 
The Hamiltonian can then be written 

L L-l 

i=l 1=1 
+ J L (c{ - c L )(c\ + ci) exp(iTrAA) , (35) 
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where 



c\ci , 



(36) 



is the number of fermions. The last term in Eq. ( |35|) is 
different from the other terms involving the Ji since the 
string operator in Eq. (|4|) acts all the way round the lat- 
tice because of periodic boundary conditions. Although 
the number of fermions is not conserved, the parity of 
that number is conserved, so exp(i7rA r ) is a constant of 
the motion and has the value 1 or —1. Hence, the fermion 
problem must have antiperiodic boundary conditions if 
there is an even number of fermions and periodic bound- 
ary conditions if there is an odd number of fermions. 
Note that the fermion Hamiltonian, Eq (|35|), is bi- linear 
in fermion operators, and so .describes free fermions. 

For the non-random modetrtl one solves for the single 
particle eigenstates of Eq. (|35| ) by (i) a Fourier trans- 
form to operators c\ and Ck, where k is the wavevector, 
followed by (ii) a Bogoliubov-Valatin transformation in 
which new fermion creation operators, 7L are formed as 
a linear combination of cl and in order to remove 
the terms in H. which do not conserve particle number. 
In the random case, we proceed in an analogous way. We 
define a column vector, 'J, and its hermitian conjugate 
row vector each of length 2L, by 



,ci,ci,c 2 , 



,c L ) 



(37) 



Note that the \I/ and W satisfy the fermion commutation 
relations 



*i*3 + *i*< 



Ijjrtlj/t 



X&tllrt 

3 1 



0, 



(38) 



irrespective of whether vPj refers to a creation or annihi- 
lation operator. For reasons that will become clear be- 
low, the Hamiltonian is written in a symmetrical form, 
replacing c;c l+ i by (cfCf+i - c i+ ic % )/2, and c\c l+1 by 
(ctcj-|-i — Cj-|_icj)/2 etc. It can then be written in terms 
of a real-symmetric 2L x 2L matrix, H, as 



U = 

where H has the form 
H = 



A B 

-B -A 



(39) 



(40) 



where A and B are L x L matrices with elements given, 
for periodic boundary conditions, by 





= hi 




= -Jill 




= -J/2 




= J/2 




= -Ji/2 , 



where i + 1 is replaced by 1 for i = L. Note that A is 
symmetric and B is antisymmetric so H is indeed sym- 
metric as claimed. For antiperiodic boundary conditions, 
one changes the sign of the terms connecting sites L and 
1 in Eq. Q. 

Next jKe diagonalize A numerically, using standard 
routineslia, to find the single particle eigenstates with 
eigenvalues e^, fi = 1, 2, . . .2L and eigenvectors $t which 

are linear combinations of the \E"J with real coefficients. 
We require that the $t have the same commutation rela- 
tions as the ^j, see Eq. (|3^), which is satisfied provided 
the transformation from the to the $ M is orthogonal, 
which in turn is guaranteed by the symmetry of H that 
we enforced above. If we interchange the c\ with the Cj 
in Eq. (|37|) then H changes sign. Hence the eigenstates 
come in pairs, with eigenvectors that are Hermitian con- 
jugates of each other and eigenvalues which are equal in 
magnitude and opposite in sign. We can therefore define 
<&t = 7^ if e M > and $L = 7^ if is the state with 
energy — e M . The Hamiltonian can then be written just 
in terms of L (rather than 2L) modes as 



n 



L 



(42) 



where all the e„, are now taken to be positive. From 
Eqs. ( |35|) and (^2|), one sees that if all the Ji are zero, 
then the e p equal the hi, as expected. We shall denote by 
"quasi-particles" excitations created by the 7^, whereas 

excitations created by the cj will be called "bare parti- 
cles". 

The many-particle states are obtained by either having 
or not having a quasi-particle in each of the eigenstates. 
One has to be careful, though, because, with periodic 
boundary conditions, the number of bare particles, j\f in 
Eq. (|3^) , must be odd, while for states with anti-periodic 
boundary conditions the number must be even. Thus, 
to generate all the many-body states one needs to solve 
the fermion problem for both periodic and anti-periodic 
boundary conditions, and keep only half the states in 
each case. 

In order to determine which states correspond to the 
the ground state and the first excited state it is useful to 
consider first the non-random caseQ'EI. There the ground 
state is in the sector with antiperiodic boundary condi- 
tions, and has no quasiparticles, which corresponds to j\f 
even as required. Hence the ground state energy is given 
by 



E 



£' 



(43) 



(41) 



where we indicate that the energies are to be evaluated 
with antiferromagnetic boundary conditions. The first 
excited state is in the the sector with periodic boundary 
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conditions. In the disordered phase, there is one quasi- 
particle, in the eigenstate with lowest energy, and this 
state has an odd-number of bare particles, as required. 
Hence the energy of the first excited state of the pure 
system in the disordered phase is given by 



Ei 



where we have ordered the energies such that t\ is the 
smallest. At the critical point of the non-random model, 
ei becomes zero. In the conventional point of view, one 
then says that t\ becomes negative in the ordered phase. 
From our perspective of numerical calculations, it is more 
convenient to define all the e M to be positive, which means 
that we are effectively interchanging the role of the cre- 
ation and annihilation operators, 7} and 71. Hence, in 
our point of view, there are now no quasi-particles, but 
this still corresponds to an odd number of bare particles. 
From either point of view, the energy of the first excited 
state of the pure system in the ordered phase is given by 



L 



(<5<0) 



(45) 



with all the taken to be positive. Note that in the 
disordered phase there is a finite gap, 2ei, in the ther- 
modynamic limit, whereas in the ordered phase the gap 
tends exponentially to zero as L — > 00. This is the man- 
ifestation of broken symmetry. Note also that we can 
rephrase the result for E\ of the pure system by saying 
that it is given by Eq. (Q) if the state with no quasi- 
particles has an even number of bare particles and by 
Eq. ([45]) , if it has an odd number (taking all the to be 
positive) . 

For the random problem the picture turns out to be 
very similar. We find that the ground state energy is 
given by Eq. ( [43| ) and the lowest excited state has en- 
ergy given either by Eq. (^4|) or Eq. (f45|), depending on 
whether the state with no quasi-jiarticles has an even or 
an odd number of bare particlealj, M . The parity of M 
is determined from Eq. (^7j) below. 

We have checked that our the code is correct by com- 
paring results for Eq and E\ for small sizes obtained from 
this fermion method with results obtained for the orig- 
inal problem, Eq. (Q), using both complete diagonaliza- 
tion and also the Lanczos method. In all cases the results 
agreed to within machine precision. 

We now proceed to the calculation of the correlation 
functions in the ground stated. As discussed above, this 
is in the sector with anti-periodic boundary conditions, 
which will be assumed in the rest of this section, unless 
otherwise stated. Assuming, without loss of generality, 
that j > i, Cij can be expressed in terms of fermions by 



Gij = (( C I +c,)exp 



i=i 



(cj + c s )) (46) 



where the averages are to be evaluated in the ground 
state. Now 



exp 



Cl 



-(4 - c i)(4 + cj) 

(4+ci){c\ -Cl) , 



(44) so defining 



Ai = c) + Q 



Bi = 



1 - c i 



and noting that A? = 1, one has 

Cij — {Bi (Ai+iBi+i . . . Aj-iBj-i) Aj) 



(47) 
(48) 



(49) 



(50) 



This rather complicated looking expression can be eval- 
uated using Wick's theorem. To see this, note first that 



(AiAj) = (6ij - c\a + c\ Cj ) 



(51) 



(since CjCi and c\cj are Hermitian conjugates of each 
other and a real diagonal matrix element is being evalu- 
ated) and similarly 



(BiBj) 



(52) 



Hence the only non-zero contractions are (AjBi) and 
(BiAj), since (BiBi) and (AiAi) never occur. Defining 



(BiAj) 



-(AjBi) = Gij , (53) 



the correlation function is given by a determinant 

Gi,i+1 Gi.i+2 ■ ■ ■ Gi 



Gi-\ 



1 G 



' ■ Gi+l J 



i+2 



•• G 



3 -hi 



(54) 



which is of size j — i. 

Gij can be expressed in terms of the eigenvectors of 
the matrix H in Eq. (Eol). Let us write 



(55) 



where ip and <p can be shown to be orthogonal matrices. 
It follows that 

G lJ = ((c!-c l )( C ]+ Cj )> 

= ^<0w((7i - 7i)(7] +7j)) 



(56) 
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since (7^7^) = (77) = and there are no quasi-particles 
in the ground state so (7^7) = 0. 

Numerically it is straightforward to compute the Gij 
from Eq. ( j56l) and then insert the results into Eq. ( 54 ) to 



determine the Cy for all i and j. 

Finally we note that the parity of the number of bare 
particles, Af, in the state with no quasi-particles can also 
be obtained, for either boundary condition, from the Gij 
since 

L 

(exp(^AO) = ([] B * A i) 

i=l 

= det G , (57) 

where we assumed that L is even, otherwise, from 
Eq. ([47]), there would be an additional minus sign. 



IV. RESULTS FOR THE ENERGY GAP 



For the pure system, the energy gap, 
AE = E 1 - E , 



(58) 



is finite in the disordered phase, and tends to zero expo- 
nentially with the size of the system in the ordered phase. 
Consider now the random case in the disordered phase, 
so [ln/i] av > [In J]av Because of statistical fluctuations, 
there are finite regions which are locally ordered, i.e. if 
one were to average just over one such region then the 
inequality would be the other way round. These regions 
will have a very small gap. Hence one expects large sam- 
ple to sample fluctuations in the gap, especially for big 
systems. 

Data for the distribution of In AE at the critical point, 
ho = 1, is shown in Fig. [l] for sizes between 16 and 128. 
One sees that the distribution gets broader with increas- 
ing system size. This is clear evidence that z — 00 as 
predicted. The precise prediction is that the log of the 
characteristic energy scale should vary as the square root 
of the length scale. With this in mind, Fig. H shows 
a scaling plot for the distribution of In AE/L 1 / , which 
works quite well. 

In the disordered phase, the data looks rather different. 
Fig. plshows the distribution of In AE, for ho = 3. Unlike 
Fig. nl the curves for different sizes now look very simi- 
lar but shifted horizontally relative to each other. This 
implies that the data scales with a finite value of z , as 
predicted. 
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FIG. 1. A plot of the distribution of the log of the energy 
gap, AE, at the critical point, ho = 1, for lattice sizes between 
16 and 128. The distribution was obtained from the value 
of the gap for 50000 samples for each size. For the larger 
sizes, the distribution is cut off at small values because, in 
this region, the gap is essentially zero within double precision 
accuracy. One sees that the distribution gets broader and 
broader as L increases. 
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(In AE) / L 1/2 

FIG. 2. A scaling plot of the data in Fig. [j], assuming that 
the log of the energy scale (here the gap) varies as the square 
root of the corresponding length scale (here the size of the 
system) . 
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FIG. 3. A plot of the distribution of the log of the energy 
gap, AE, in the disordered phase at ho = 3, for lattice sizes 
between 16 and 128. The distribution was obtained from the 
value of the gap for 50000 samples for each size. The different 
curves are very similar and just shifted relative to each other. 

Note that in the region of small gaps, the data in Fig. || 
is a straight line indicating a power law distribution of 
gaps. This power law behavior is not special to the l*d 
problem discussed here, but is expected quite generallyO 
in the Griffiths phase for systems with discrete symmetry. 
The power is related to z as we shall now see. Well into 
the disordered phase, excitations which give a small gap 
are well localized so we assume that the probability of 
having small gap is proportional to the size of the system, 
L. This assumption is confirmed by the data in Fig ^. 
Hence, the probability of having a gap between AE and 
AE(1 + e) (for some e) should have the scaling form, 
eLAE 1 ^, so the distribution of gaps, P(AE), must vary 
as 

P(AE) ~ AE~ l+1 ' z , (59) 

in the region of small gaps. It is tidier to use logarith- 
mic variables, and the corresponding expression for the 
distribution of In AE is 

In [P(ln AE)] = -\nAE + const. (60) 

z 

From the slopes in Fig. || we estimate z ~ 1.4,which 
gives a satisfactory scaling plot as shown in Fig. [|. The 
data does not collapse so well for large gaps, but this may 
be outside the scaling region. 
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In (L z AE) 

FIG. 4. A scaling plot of the data in Fig. |^, assuming scal- 
ing with a finite value of z. The fit here has z = 1.4. 




0.5 1 

6 

FIG. 5. Results for 1/z against 5, where z is the dynamical 
exponent and 8 is related to ho by Eq. ([]) . The solid curve is 
a fit fcp 1/z = 28( 1 — 2SC) (which corresponds to the expected 
formQ, Eq. (|)) with C = 0.311. The dashed line is l/z = 25, 
the predicted asymptotic form for 5 — + 0. 
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FIG. 6. A plot of the distribution of the log of the energy 
gap, AE, at the critical point, ho = 1, for the bimodal distri- 
bution in Eq. (plj). The distribution was obtained from the 
value of the gap for 50000 samples for each size. 

We have carried out a similar analysis for other values 
of ho. Close to the critical point, it is difficult to deter- 
mine z because the distribution broadens with increas- 
ing size for small sizes (presumably where L < £), but 
then the slope of the straight line region starts to satu- 
rate, corresponding to a large but finite z. The sizes that 
we can study are therefore in a crossover region between 
conventional dynamical scaling (z finite) and activated 
dynamical scaling [z infinite) so the data does not scale 
well with, any choice of z. 

Fishem has predicted that z is equal to 1/2S + C, near 
the critical point, where C is non-universal constant, see 
Eq. (||) . We show our estimates for l/z plotted against 6 
(which is related to ho by Eq. (||)), in Fig |[ Also shown 
is a fit of l/z to 25(1 — 25C) , which corresponds to Eq. (||) 
and which works quite well with C = 0.311. 

The exponents are predicted to be universal, i.e. in- 
dependent of the distributions 7r(J) and p(h) (as long as 
these don't have anomalously long tails). To test uni- 
versality, we also did some calculations for a bimodal 
distribution, in which J and h take one of two values, 



n(J) = ^[6(J-l) + 6(J-3)] 
p(h) = ~[5(h-ho)+8(h-3h )] 



(61) 



The critical point is at ho = 1. The data for AE at the 
critical point is shown in Fig. |^ and the scaling plot is 
presented in Fig. 0. The data scales reasonably well in- 
dicating that z — oo at the critical point, just as for the 
continuous distribution in Eq. (|^). The data collapse is 



not as good as for the continuous distribution, presum- 
ably indicating that the approach to the scaling limit is 
slower. 
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FIG. 7. A scaling plot of the data in Fig. |^. The collapse of 
the data indicates that z — oo, the same as for the continuous 
distribution of Eq. see Fig. ^. 



V. RESULTS FOR CORRELATION FUNCTIONS 

We start by looking at the correlation functions at the 
critical point and then discuss our results in the disor- 
dered phase. 

The average correlation function at the critical point is 
shown in a log- log plot in Fig || for several sizes. The data 
for the larger sizes lie on a straight line, and the dashed 
line, which is a fit to the L = 128 data for 7 < r < 35, has. 
a slope of —0.38, in excellent agreement with Fisher'scl 
prediction in Eq. (|l4|). 

A graph of the average of the log of the correlation 
function (which corresponds to the log of a typical corre- 
lation function) is shown in Fig. ^| plotted against y/r. As 
expectedEm from Eq. (22) the data falls on a straight line. 
The data in Figs. || and ^| indicate that the average and 
typical correlation functions do behave very differently 
at the critical point, as predicted. 
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FIG. 8. A log-log plot of the average correlation function 
against distance at the critical point. The straight line be- 
havior for larger sizes indicates a power law variation. The 
dashed line is a fit to the data for L — 128 with 7 < r < 35 
and has slope of —0.38 in excellent agreement with the pre- 
diction in Eq. (|l4|). The results are obtained by averaging 
over all pairs of points separated by a distance r for 10000 
samples. 




FIG. 9. A plot of the average of the log of the correlation 
function against \fr at the critical point. Thiftrstraight line 
behavior for larger sizes supports the predictionEm of Eq. (p2[). 
The results are obtained by averaging over 10000 samples. 
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FIG. 10. The distribution of the log of the correlation 
function for different values of r at the critical point. The 
data is obtained from 10000 samples of size L — 128. 
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FIG. 11. A scaling plotrtpf the data in Fig. |l0| according 
to the theoretical predictions. The data collapses well in the 
region of fairly high probability, including the upturn near 
the right hand edge (which may indicate a divergence as the 
abscissa tends to 0, see the text). There are systematic devi- 
ations in the tail which may be due to corrections to scaling. 

The reason for this difference is that the distribution 
of In C(r) is vepy broad, as can be seen in the plot in 
Fig. [l(]. Fishem has predicted that the distribution of 
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(In C(r))/y/r should be universal and independent of r so 
we show the corresponding scaling plot in Fig. O. Note 
that the distribution monotonically decreases as C(r) be- 
comes smaller. The data scales well for larger values of C, 
including even the upturn near the right hand edge of the 
graph. This is the region where C(r) is anomalously large 
and which gives the dominant contribution to the aver- 
age correlation function. An interesting question, then, 
is whether the value for the average correlation function 
is included in the scaling function for In C(r)/y / r. If so, 
the the^scaling function will diverge as a power near the 
origin!!! To see this, note that if the probability of hav- 
ing a correlation C at a distance r only depends on the 
combination 



and that if 



j/=(lnC)/V^ 



P(y) 



(62) 



(63) 



An enlarged log-log plot of the region of the upturn 
in Fig. [n]. is shown in Fig. [l^. The data does lie on a 
rough straight line whose slope decreases (in magnitude) 
with increasing r. A fit to the data for r = 24 in the 
middle region of the graph has a slope of about 0.45 (in 
magnitude), larger than 0.24, but since the effective slope 
decreases with increasing r the data does not rule out the 
possibility that the distribution of (In C(r))/y/r diverges 
with an exponent of 0.24 for r — > oo. Even if the scaling 
function for (In C(r))/y/r gives the correct power for the 
average correlation function, it does not necessarily mean 
that the amplitude is correct, since there could also be 
additional non-universal contributions to the amplitude, 
outside the scaling functionlill. 

We suspect that the systematic deviation in the tail 
of the distribution in Fig. |ll| at small values of C(r), 
indicates corrections to scaling for this range of sizes and 
distances. 



for small y, then the average value of C(r) is given by 



[C(r)} £ 



dC 



lnCV A 1 



r-l/2 



„l/2 ' 



(64) 



assuming that the integral is dominated by the region of 
small In C . Integrating over C gives a finite number so 



and comparing with Eq. ( |i"4| ) yields 
A = 20 - 3 ~ 0.24 . 



(65) 



(66) 
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-(In C(r)) / r 1 / 2 

FIG. 12. An enlarged log-log plot of the region of the 
upturn in Fig. [ll] where the abscissa approaches zero.. 
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FIG. 13. A scaling plot of the average correlation function 
in the disordered phase according to Eq. (po[). The^best fit 
(shown) is for v — 1.8, fairly close to the predictions v — 2. 
The data represents an average over 10000 samples. 
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FIG. 14. A sca ling plot of the average correlation function, 
according to Eq. (fecj) with the predicted value v = 2. The 
curve is a plot of Eq. ( [ill ) (expected to be valid for rti 2 2> 1) 
with D = 35. Note that from Eqs. (O) and (Jl9|) , one has 
= S 2 for the distribution of Eq. 
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FIG. 15. A scaling plot of the average of the log of 
the correlation function in the disordered phase according to 
Eq. (EfiHnThe best fit (shown) is for v = 1.1, close to the pre- 
dictionQ'S v = 1. The data is an average over 10000 samples. 





-1 

■2 
3 
■4 
-5 
■6 
-7 



iiiiiiiii nun 


II II III 1 II 1 II l : III 




"iirrff 


= \, 


h 




- 


- 


1.3 


• 


- 




1.5 


X 




= z 


it 1-7 


X 




= 


r£ 2.0 


«t 


- 




^% 2.3 


O 


= 


= 














= 






































\x 












1 


ii 1 1 nil 


nil ik 


in i ft 



12 3 



4 5 6 7 
2r6 



FIG. 16. A scaling plot of the average of the log of 
the correlation function in the disordered phase according to 
Eq. (p^), assuming the theoretical value = 1. The solid line 
is the prediction of Eq. (pTj), which expected to be valid at 
large r, but also works well down to r = 0. It appears that 
by dividing by the correlation function at criticality, we in- 
corporate most of the corrections to the asymptotic form in 
Eq. @. Note that from Eq. @ one has f" 1 = 25. 

We now discuss our results for the disordered phase. 
The scaling plot corresponding to Eq. ( |20| ) is shown in 
Fig. [l3|. The plot has v — 1.8 which gave the best fit, 
and which agrees fairly well with the prediction v — 2. A 
plot using v = 2 works somewhat less well, presumably 
indicating that there are corrections to scaling for this 
range of lattice sizes and distances. Fig. [l4] is a scaling 
plot using the theoretical value v = 2 which also shows 
the asymptotic form in Eq. ( pl| ) with D = 35. Both 
the data and the prediction of Eq. (|2l]) , have substantial 
curvature: much more than in the corresponding data 
for log C(r) shown in Figs. |l5| and [l(| Over the range 
of accessible values of rS 2 , the data and the asymptotic 
prediction do not track each other closely though it is 
possible that they would do so for larger values of r5 2 . 

The data for the log of the scaling function scales well 
according to Eq. ( |26| ) , though the best fit has a slightly 
different exponent of 1.1, see Fig. |l5|. Presumably this 
difference again indicates that there are corrections to 
scaling for the sizes and distances studied. Note that the 
data in Fig. [15] is close to a straight line but there is 
statistically significant curvature. |-.|— . 

Fig. |l6| tests the more stringent predictionBO for the 
average of In C (r) in the limit r £ obtained by com- 
bining Eq. (gTj) with the assumption that the expression 
for | in Eq. (ps|) is exact for correlations of a z . One sees 
that it works quite well. 
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FIG. 17. The distribution of the log of the correlation 
function for different values of r at ho = 3 in the disordered 
phase. The data represents an average over 10000 samples 
for L — 64. Both the most probable value and the width of 
the distribution increase with r but the most probable value 
increases faster, see Fig. [j^. 
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FIG. 18. A scaling plot of the data in Fig. |T?]. The width 
of the distribution is seen to narrow with increasing r, as 
expected. 

Although our best estimates of the critical exponents v 
and v do not quite agree with the theoretical predictions, 
they are fairly close to those predictions, and they differ 



substantially from each other, providing clear evidence 
that there are different correlation length exponents for 
the average and typical correlation functions. 

Finally, in this section, we look at the distribution of 
In C(r) for r larger than either tha. average or typical 
correlation lengths. It is predictedEm that the distribu- 
tion of (In C(r))/r should become sharp at large r in this 
limit. Fig. 17 shows data for the distribution of In C{r) 
at ho — 3. One sees that both the peak position and the 
width increase with increasing r, but the peak position 
increases faster as can be seen in Fig. |l8| w hich shows 
the distribution of (In CLf^/r. In Fig. |19| we test the 
more precise predictionsuEil for the mean and variance 
of lnC(r) given in Eqs. @) and @, for ho = 3. The 
fits give reasonable agreement with Eqs. (^lj) and (|3^), 
but assume a form of corrections to scaling that we have 
been unable to justify. 
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FIG. 19. The mean and variance of the distribution of 
In C(r) in the disordered phase at ho — 3 for different values of 
r. Both the mean and variance are expected be proportional 
to r at large r. The lines are least squares fits of the form 
a + br 1 / 2 + cr. This form is motivated by the behavior at 
criticality, for which the mean varies as r 1 ^ 2 , and, in fact, the 
r 1//2 correction to scaling was effectively removed in Fig. |l6| 
by factoring out the behavior at criticality. However it is 
not clear that the constant and the r 1 ^ 2 term give all the 
corrections to scaling in the disordered phase. The fit to the 
mean has a = —0.140, b = 1.145, c = 1.014, while the fit to the 
variance has a = 1.025, fe-p — 1.736, c = 2.029. According to 
the suggestion of FisheiO, following Shankar and Murphyu, 
the leading behavior for large r should be given by Eqs. ( |3l| ) 
and (^). Noting that here, 25 — In 3 ~ 1.099, we see that 
the large r behavior predicted by the fits is is in rather good 
agreement with theory. 
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VI. THE LOCAL SUSCEPTIBILITY 

In this section we discuss the local susceptibility rather 
than the uniform susceptibility, because it has somewhat 
simpler behavior. Since it it just involves correlations 
on a single site, any singularity must come only from 
long time-correlations, whereas the uniform susceptibil- 
ity involves correlations both in space and time. Our 
results for the uniform susceptibility away from the criti- 
cal point do not scale in a simple manner, and we suspect 
that there are logarithmic corrections, as occurs for bulk 
behavior at finite temperatureQ. 
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FIG. 20. The distribution of the log of the local suscepti- 
bility at the critical point for different sizes obtained by the 
Lanczos method. The number of samples is 50000/L, so that 
50000 values for Xioc were obtained for each size. 

Since it is difficult to compute the susceptibility from 
the fermion method, particularly with periodic bound- 
ary conditions, we have used the Lanczos diagonaliza- 
tion technique on the original spin Hamiltonian, Eq. (Q). 
Of course the price we pay is that the lattices are much 
smaller, L < 16. The local susceptibility at T = is 
given by 



Xlo 



E n — Eq 



(67) 



where \n) denotes a many body state of the system and 
|0) is the ground state. Because of the form of Eq. j67| ) 
we expect that the scaling of % loc will be very similar to 
that of l/AE. This is indeed the case as seen in Fig. pc| , 
which plots the distribution of lnx loc at the critical point. 
The distribution broadens with system size, consistent 
with z — oo. The data scales in the expected manner, 



as shown in Fig. gl], which is very similar to the corre- 
sponding plot for the energy gap in Fig. ||. 
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FIG. 21. A scaling plot of the data in Fig. E3. The data 
scales in the same way as that for the log of the gap, see Fig. |j. 
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FIG. 22. The distribution of the log of the local suscep- 
tibility in the disordered phase for different sizes obtained by 
the Lanczos method. The number of samples is 50000 / L. 

Even though the range of sizes used in the Lanczos 
method is rather small, it is, nonetheless, capable of dis- 
tinguishing z — oo scaling at the critical point from finite 
z scaling z away from the critical point. This can be seen 
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by comparing Fig. |2(J with Fig. g2], which plots the dis- 
tribution of lnx loc at ho = 3. In Fig. 22 the curves no 
longer broaden with increasing L but the distributions 
are independent of size. The reason why there is no size 
dependence here but there is in the distributions of In AE 
in Fig. H is easy to understand. For AE, we compute the 
probability per sample of getting a certain value, and this 
is proportional to L in the disordered phase for small AE, 
since the rare strongly correlated region can occur any- 
where. We used this result in Section IV to relate the 
exponent in the distribution to 1/z, see Eqs. ( p9[ ) and 
(|60|). With Xioc however, we compute the probability 
per site, so there is no factor of L and the distribution is 
independent of size. This is, of course, the normal state 
of affairs when the lattice size is much larger than the 
correlation length. The slope of the straight line region 
in Fig. 22 agrees with the slopes in Fig. || and so gives the 
same value of z as obtained from the gap, i.e. z ~ 1.4. 



VII. THE STRUCTURE FACTOR 

A scattering experiment measures directly the struc- 
ture factor, S(q), defined by 



(68) 



3,1 



Although the distribution of individual terms in the 
sum is broad, it is interesting, and relevant for experi- 
ment, to ask whether there are large sample to sample 
fluctuations in the total. We have attempted to answer 
this for q = 0, where fluctuations are expected to be 
largest. Fig. ^3| shows a log-log plot of the average of 
the structure factor, S av (0), and the standard deviation 
among different samples, (55(0), plotted against L at the 
critical point. From Eq. (|lj) one expects the average 
to vary as L 62 and the best fit to the numerical data 
has a slope of 0.64, in reasonably good agreement. One 
sees that 5S(0) < S av (0), but the c&jtio of the width to 
the mean stays finite. One expectsM that the distribu- 
tion of 5(0)/5 av (0) will be broad and independent of L 
at large L, and our data is consistent with this. Note 
that although the equal time structure factor is not self- 
averaging at T = 0, its distribution is much less broad 
that that of the susceptibility, which involves correlations 
in time. Since the structure factor at the critical point 
behaves like the average, rather than the typical corre- 
lation function, we expect that the behavior away from 
criticality will be controlled by the average correlation 
length. 
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FIG. 23. A log-log plot of the average q = structure 
factor, 5 , av (0) and the standard deviation of the structure 
factor among different samples, ^(O), for different sizes. The 
dashed line is a fit to the data for 5' av (0) and has slope 0.64, 
in fairly good agreement with the theoretical expectation, ob- 
tained by integrating Eq. (|l4|), of (f>— 1 ~ 0.32. The solid line 
is just a guide to the eye. The results are obtained by aver- 
aging over 10000 samples. 



VIII. CONCLUSIONS 

We have been able to confirm the many surprising pre- 
dictions of the random transverse-Ising spin chain by ap- 
plying the mapping to free fermions numerically. In par- 
ticular we find very broad distributions of the energy gap 
and correlation functions, different exponents for the av- 
erage and the typical correlation functions, and an infi- 
nite value of the dynamical exponent, z, at the critical 
point. Perhaps the most interesting new result is the 
scaling function for the distribution of the logof the cor- 
relation function at criticality, shown in Fig. |ll| , which is 
monotonic and has an upturn as the abscissa approaches 
zero. If this indicates the divergence shown in Eqs. ( |6^ ) 
and (|66|), the scaling function for (In C(r))/r 1 / 2 would 
also give the correct exponent (though perhaps not the 
correct amplitude) for the average correlation function. 

We have seen that the width of the distribution of the 
the equal time structure factor seems to be comparable 
to the mean at T = 0, though presumably it becomes 
self-averaging at finite-T. By contrast, the T = suscep- 
tibility and local susceptibility have enormously broad 
distributions. One expects that the susceptibility will 
also become self averaging at finite T for sufficiently large 
L, but whether the necessary size diverges as power law 
or exponentially as T — ► is unclear. We leave this in- 
teresting question for future study. 
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Crisanti and RiegerB have studied the random trans- 
verse Ising chain by Monte Carlo methods. They took 
a generalization of the bimodal distribution in Eq. j6l| ) 
rather than the continuous distribution used here. From 
the behavior of various correlation functions they found 
a finite z at criticality, which, however, appeared to in- 
crease with increasing randomness. We saw in section IV 
that corrections to finite-size scaling appear to be larger 
for this distribution than for the the continuous one and, 
furthermore, it is harder to estimate the asymptotic value 
of z from correlation functions than from distributions. 
This is presumably why Crisanti and Riegerta did not 
find z = oo in their study. 

After this work was largely completed we .became 
aware of related work by Asakawa and Suzuki£L3, who 
also used the mapping to free fermions but used the same 
distribution as Crisanti and Rieger. In contrast to our 
results, they claim that the exponents depend on the pa- 
rameters in the distribution- This is lack of universality 
is not predicted by theorytffl and a possible explanation 
of the discrepancy is that not all their data is in the 
asymptotic scaling regime, which is likely to be reached 
for different lattice sizes for different distributions. 

It is interesting to speculate to what extent the re- 
sults of the one-dimensional system go over to higher 
dimensions. In particular, one would like to know if z 
is infinite at the critical point or takes a finite value for 
d > 1. The results for the local susceptibility in Section 
VI indicate that this question can be answered even for 
moderately small lattice sizes provided appropriate quan- 
tities are studied. The distribution of In x loc (or the log 
of the gap) seems to be particularly convenient, since, for 
finite- z, data for different sizes look essentially the same, 
whereas for z = oo the curves get broader and broader. 
Of course it is still difficult to distinguish a large but fi- 
nite z from z — oo, since the two would look the same 
for small sizes. With finite-z scaling, the distribution 
has power law behavior, the power being related to z 
as shown in Eq. (60). It is more difficult to determine 
z by looking at the decay of correlation functions, be- 
cause the asymptotic behavior is only seen at very large 
times or distances. Numerical studies in higher dimen- 
sions are likely to use quantum Monte Carlo simulations, 
because diagonalization methods, such as Lanczos, can 
only be carried out on very small systems and the map- 
ping to free fermions only works in one-dimension. Un- 
fortunately, there is an additional difficulty with quan- 
tum Monte Carlo, not present here, because one gener- 
ally works in imaginary time, which has to be discretized. 
The quantum problem is recovered when the number of 
time slices tends to infinity, but in practice one can only 
simulate a finite number. It is unclear whether the ex- 
trapolation to an infinite number of time slices will pose 
serious difficulties for the study of Griffiths singularities 
and critical phenomena in higher dimensional systems. 

We have seen that the disordered Griffiths phase can 
be conveniently parameterized by a continuously varying 
dynamical exponent z. This characterizes the distribu- 



tion of the energy gap or local susceptibility for lattice 
sizes which satisfy the condition, L £. By contrast, at 
the critical point, the correlation length diverges so the 
value of z at criticality involves physics in the opposite 
limit, L <C £. It is therefore possible that the limit of 
z{8) for 8 — > is not equal to the value of z at criticality. 
Both these quantities are infinite for the transverse field 
Ising chain, but it would be interesting to see if there 
is a difference between them in higher dimensions. We 
expect that the results and method of analysis presented 
here will provide guidance for such a study. 
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For the random problem, as well as for the pure problem, 
one eigenvalue hits zero as ho (or equivalently S) is varied. 
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tered around the critical point. Hence, for the pure system, 
one knows to choose Eq. ( fli] ) or ( p5[ ) for E\ depending on 
whether the system is in the disordered or ordered phase, 
but there is no such simple criterion for the random case. 
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